function [ params ] = cellParamsFromVsteps( varargin )
%CELLPARAMSFROMVSTEPS Calculates cell parameters from current response to
%voltage steps (in voltage-clamp mode)
%
% Uses the method described in the pClamp manual to calculate Rm, Rs, Cm,
% tau and I_hold
%
% % Optional Arguments
%  1: showResponseFig boolean: show/hide the averaged response figure
%  2: currentTrace struct (if not present, will load from gcf.UserData)
%
% Ben Suter, 2009-09-26
% Ben Suter 2009-09-12 -- added output and input args, conditional figure
% display, proper output format, etc. to match the latest "standard"
% version; also prints results ready for copy-paste into .m file
% Ben Suter 2009-09-12 -- fixed bug in output: I_hold sign was inverted
% Ben Suter 2009-11-09 -- fixed output bug, was printing bl not I_hold
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If no currentTrace argument, load from UserData of current figure
if ( nargin < 1 )
    showResponseFig = true;
    data = get(gcf, 'UserData');
    currentTrace = data.currentTrace;
else
    showResponseFig = varargin{1};
    currentTrace = varargin{2};
end

stimulus = currentTrace.header.ephys.ephys.pulseParameters{1};
if isempty(stimulus)
   error('Unable to get pulse parameters (might not be item 1?)');
elseif ~ strcmp(stimulus.type, 'squarePulseTrain')
   error('Stimulus must be a square pulse train');
end

%                       type: 'squarePulseTrain'
%                 sampleRate: 10000
%                  amplitude: -5
%                     offset: 0
%     squarePulseTrainNumber: 4
%        squarePulseTrainISI: 0.2500
%      squarePulseTrainWidth: 0.1000
%      squarePulseTrainDelay: 0.1000

delay = stimulus.squarePulseTrainDelay;
ISI = stimulus.squarePulseTrainISI;
width = stimulus.squarePulseTrainWidth;
numPulses = stimulus.squarePulseTrainNumber;
amplitude = stimulus.amplitude;
sampleRate = stimulus.sampleRate;

if delay < width / 2
   % We use 1/2 the pulse width preceding the pulse as baseline
   % and the 2nd half of the pulse width for the steady-state
   error('The delay before the first pulse must be no less than half the pulse width');
end

trace = currentTrace.data.ephys.trace_1;

% Snippet: 1/2 pulse width preceding stimulus onset (for baseline)
% and 1 pulse width post stimulus onset (for response)
for i=1:numPulses
   % only using negative-going response for now
   % TODO: look at return-to-baseline response as well
   rs = 1 + (delay + (i-1)*ISI)*sampleRate; % response start
   re = rs + width * sampleRate - 1; % response end
   response(i,:) = trace(rs:re);
   baseline(i,:) = trace(rs-width/2*sampleRate:rs-1);
   steadyState(i,:) = trace(re-width/2*sampleRate+1:re);
end

% average over all pulses
resp = mean(response, 1);
bl = mean(mean(baseline, 1));
ss = mean(mean(steadyState, 1));

% I_hold is baseline
I_hold = bl;

if ( amplitude < 0 )
   resp = - resp;
   bl = - bl;
   ss = - ss;
end

% t = 0:width*sampleRate - 1;
t = 0:1/sampleRate:width - 1/sampleRate; % time from zero in seconds

if showResponseFig
    figure; hold on;
    plot(t, resp);
    plot(0, bl, 'g*'); plot(0, ss, 'r*');
    plot(t, resp-ss, 'y');
end

% tau:
% Fit exponential decay to response, between 10% and 80% of the
% peak response (these numbers are taken to match the pClamp "Proportion
% of peak to fit" defaults)
respExp = resp - ss; % response above steady-state, i.e. exponential decay
[peak peakIndex] = max(respExp);
% Sometimes the first sample post-stimulus is only part-way up the rise,
% i.e. not the peak value. I assume this is due to the way ephus or the
% DAQ boards digitize the signal (average, interpolate?). The fitting will
% be easier if we truncate any values preceding the peak value.
respExp = respExp(peakIndex:end);
t = t(peakIndex:end);
% Start below 80% and stop before 10% of peak decay
fitStart = min(find(respExp < 0.8*peak));
fitStop = max(find(respExp > 0.1*peak));
% Fit a single exponential to the decay
[estimates, model] = fitExp(t(fitStart:fitStop), respExp(fitStart:fitStop));
tau = 1000 / estimates(2); % tau, converted to ms

% Q1 = Integrate I above steady-state
q1 = sum(resp - ss) / (sampleRate / 1000); % pA / samples per ms

% Q2: The amount of charge below steady-state and above the resistor
% current
q2 = tau * (ss - bl);  % ms * pA

% all in ms, mV , pA, should result in GOhm
Rs = tau * abs(amplitude) / (q1 + q2);
Rm = abs(amplitude / (ss - bl)) - Rs;
Cm = (q1 + q2 ) / abs(amplitude) * (Rm + Rs) / Rm; % in pF

Rs = Rs * 1000; Rm = Rm * 1000; % Convert to MOhm

if showResponseFig
    % Plot the region to be fit and the fitted exponential
    plot(t(fitStart:fitStop), respExp(fitStart:fitStop), 'r');
    plot(t(fitStart:fitStop), estimates(1) .* exp(-estimates(2) * t(fitStart:fitStop)), 'c');
    % xlim([0 max(t)]);
    xlim([0 2*t(fitStop)]);
end

% command line display
disp('Rs   Ri   Cm   tau   I_hold');
disp([num2str(round(Rs)) '  ' num2str(round(Rm)) '  ' num2str(round(Cm)) '  ' num2str((tau/1000)) '   ' num2str(round(I_hold))]);

% results -> workspace
cellParameters.Rs = Rs;
cellParameters.Ri = Rm;
cellParameters.Cm = Cm;
cellParameters.tau = tau / 1000;
cellParameters.I_hold = I_hold;
assignin('base', 'cellParameters', cellParameters); 

params = cellParameters;

% %%%%%%% Optionally print results in .m-file ready format %%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
doDataPrint = true;
if doDataPrint
    prefix = 'map.vsteps.';
    disp('% ********** V-steps: cell parameters **********')
    fnames = fieldnames(params);
    for n = 1 : size(fnames,1)
        disp([prefix fnames{n} ' = ' num2str(params.(fnames{n})) ';']);
    end
end
end